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ABSTRACT 

Although cosmological solutions to Einstein's equations are known to be generically 
singular, little is known about the nature of singularities in typical spacetimes. It is 
shown here how the operator splitting used in a particular symplectic numerical integration 
scheme fits naturally into the Einstein equations for a large class of cosmological models 
(whose dynamical variables are harmonic maps) and thus allows study of their approach 
to the singularity. The numerical method also naturally singles out the asymptotically 
velocity term dominated (AVTD) behavior known to be characteristic of some of these 
models, conjectured to describe others, and probably characteristic of a subclass of the 
rest. The method is first applied to the generic (unpolarized) Gowdy T"^ cosmology. Exact 
pseudo-unpolarized solutions are used as a code test and demonstrate that a 4th order 
accurate implementation of the numerical method yields acceptable agreement. For generic 
initial data, support for the conjecture that the singularity is AVTD with geodesic velocity 
(in the harmonic map target space) < 1 is found. A new phenomenon of the development 
of small scale spatial structure is also observed. Finally, it is shown that the numerical 
method straightforwardly generalizes to an arbitrary cosmological spacetime on x R 
with one spacelike U{1) symmetry. 



I. Introduction. 

Powerful theorems [1] prove singularities to be a generic feature of Einstein's equa- 
tions yet say nothing about the nature of these singularities. In particular, little is known 
about the singularity behavior of generic spatially inhomogeneous cosmologies. Belinskii, 
Khalatnikov, and Lifshitz (BKL) [2] and coworkers [3] have long argued that the Mixmas- 
ter dynamics [2, 4] of spatially homogeneous Bianchi Type VIII and IX cosmologies [5] 
characterizes the generic "big bang." Their results are not generally accepted, however [6], 
and evidence suggests that Mixmaster behavior disappears in models with more than three 
dynamical degrees of freedom [7] . An alternative to Mixmaster dynamics is asymptotically 
velocity term dominated (AVTD) behavior where (heuristically) terms in Einstein's equa- 
tions containing spatial derivatives can be neglected in favor of those with time derivatives 
[8, 9]. Near the singularity, AVTD solutions can be interpreted as a different spatially ho- 
mogeneous cosmology at each point in space. The polarized Gowdy cosmologies [10] have 
been shown rigorously to belong to this class [9, 11, 12]. It has recently been conjectured 
that the general (unpolarized) Gowdy models are AVTD [13]. 

We propose to study the approach to the singularity numerically using a method 
uniquely suited to the task. For both the Gowdy cosmologies defined to have a symmetry 
plane (i.e. two spatial, hypersurface-orthogonal, surface-forming Killing fields) and the 
more general cosmology possessing a single, spatial U{1) symmetry [14], (at least some of 
the) degrees of freedom can be understood as harmonic maps [15]. The superhamiltonian 
whose variation yields these equations is just an energy- like expression of the harmonic 
map fields [16]. The variation of the "kinetic" term alone yields the AVTD equations of 
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motion for these fields (with all spatial-derivative-containing terms obtained upon variation 
of the "potential" term) . This suggests that a symplectic scheme for numerical integration 
that separately evolves the kinetic and potential energy operators to approximate the total 
hamiltonian evolution is ideally suited to this problem [17]. In the following discussion, we 
shall demonstrate that this is indeed the case. 

The approach to the singularity of the Gowdy T^ cosmology has been used to test 
the feasibility of our approach. In the process of code development and testing, we have 
been able to demonstrate AVTD behavior of the singularity approach for several sets of 
(presumably) generic initial data. We have also noted that the solutions develop a char- 
acteristic small scale spatial structure which represents a competition between nonlinear 
growth and the approach to the AVTD regime which freezes the spatial profile of the 
wave amplitudes. This richness of the Gowdy T^ phenomenology will be discussed else- 
where [18]. Here we wish to emphasize the applicability of our methods to the study of 
cosmological singularities. 

While the numerical study of plane symmetric cosmologies has been underway since 
the late 1970's [19], previous work tended to focus on analogies between the cosmologi- 
cal problem and the original numerical studies of colliding black holes [20]. This led to 
concentration on constant mean curvature foliations and on the choice of lapse and shift. 
Physical interest centered on interacting wavepackets of a single polarization. In contrast, 
we begin with the predefined coordinate system in which the equations are known to be 
relatively simple and within which the harmonic map structure can be seen. This foliation 
naturally selects (in the plane symmetric case) a measure of area in the symmetry plane as 
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the time variable. Our approach then easily allows study of unpolarized Gowdy models 
and appears to be generalizable to the ^7(1) problem. 

In Section II, we shall describe the numerical method including its generalization to 
arbitrary order of accuracy in both time [21] and space. To date, our primary application 
of this method has been to the approach to the singularity of the unpolarized Gowdy 
cosmology [10, 12, 16, 22]. In Section III, the relevant properties of this model will be 
reviewed. Section IV demonstrates the validity of the code with a "pseudo-unpolarized" 
test model constructed by boosting in the harmonic map target space an analytic solution 
for the polarized model. In Section V, we discuss results for a generic unpolarized model, 
demonstrating AVTD behavior and briefly discussing the appearance of small scale spatial 
structure. In Section VI, the applicability of this method to the U{1) problem is shown. 
A summary and conclusions are given in Section VII. 
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II. The Symplectic Integrator (SI) [17] 

For convenience, we shall restrict our discussion of the method to a single degree of 
freedom which depends on only one spatial dimension and time — q{x,t) and its conjugate 
momentum 7r{x,t). We assume that the equations of motion can be obtained by variation 
of a hamiltonian 

H = ^dx[^7r^ + V{q)]. (2.1) 
Consider a differenced form of (2.1) 



N r 2 



i=0 L 



(2.2) 



where we assume periodic boundary conditions on the lattice with labels denoting 
the point {xi, P). The potential Vi at the point Xi may depend on the value of q at several 
spatial grid points. 

The symplectic scheme splits the hamiltonian operator as 

H = H^ + H2 (2.3) 



where 

Hi = jdx Itt^ (2.4a) 

and 

H2 = jdxV{q) (2.46) 

respectively. It is convenient to represent the scheme using quantum mechanical notation. 
It is based on the second order (in the time-step e) accurate approximant to the evolution 
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operator: 

^-isH ^ ^-{i/2)eH,^-^eH,^-{i/2)eH, ^ ^^^3) (2.5) 

i.e. to evolve (tt^, qi) at {xi,P) to (tt^"*""^, qj'^^) at {xi, P +e), evolve with H2 for 1/2 time- 
step, with Hi for a full time-step, and with H2 for 1/2 time-step using the appropriate 
intermediate result at each stage. In this evolution. Hi and H2 are separately to be 
regarded as the hamiltonian of the system. In the case where Hi and H2 can be separately 
exactly solved, the implementation of the method becomes trivial. 
For the hamiltonian (2.1), the scheme becomes [17] 



qi^' = qi + e 



(2.6) 



where Fi{qi^) — —dV j dgj is the appropriate force component. As an example, we consider 
the wave equation with 

where A is the lattice spacing in the x direction. Direct substitution shows that, in this 
case, the method is equivalent to the standard leap-frog differenced form of the wave 
equation [23] 

Oi^' + Oi-" - 4+1 - 4-1 = 0- (2.8) 

However, the SI algorithm has significant advantages over the leap-frog scheme for our 
problem of the approach to the singularity. 



1. As a symplectic scheme, the evolution takes the form of a canonical transformation from 
the beginning to the end of the time-step [17]. This may help to preserve the constraints of 
the cosmological problem during the evolution. (Although the continuum Einstein equa- 
tions automatically preserve the constraints during evolution, there is no corresponding 
statement for the discretized equations. Ultimately, this is a consequence of the role of 
the constraints as the generators of diffeomorphism invariance — a fundamental property 
of the continuum.) 

2. In the cosmological case. Hi is the hamiltonian whose variation yields the AVTD 
equations. If the solution is in the AVTD regime, then this SI will become increasingly 
more accurate. 

3. To avoid the problematic need to solve and resolve the constraint equations at frequent 
intervals [24, 20], one can try to find an accurate solution to the dynamical equations. 
Suzuki has shown how to generalize the SI to an arbitrary order in time [21]. The idea is 
to find an approximant like (2.5) accurate to the desired order. Such a program does not 
have a unique solution. Since it can be shown [21] that a 2m — 1 order accurate scheme is 
also 2m order accurate, one finds the recurrence relation [21] 



S2m-l{£) — S2m{^) — S2m.-z{km^)S2m-z[{^ — 2km)s\S2m-Z 



(/CjTjS) 



(2.9) 



where 




(2.10) 



with 




(2.11) 
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Thus the higher order scheme can be constructed from time-steps of the appropriate du- 
ration of the second order scheme. 

The generahzation of the spatial evolution to arbitrary order is simple only for one 
spatial dimension. For higher dimensions, the construction must proceed on a case by case 
basis. In the one spatial dimension wave equation, the second derivative with respect to x 
can be obtained in a differenced form (for spatial grid interval A) as 



The coefficients are chosen to cancel the terms in the Taylor expansion containing higher 
derivatives to the indicated order. The expression on the right hand side of (2.12) can be 
obtained as the negative of the variation with respect to fi (where fi = f{qi)) of 

N n-1 

= (2-13) 

i=0 fe=l 

where the ak are the same as those in (2.12) and with = —2 '^fe- Thus the second 

order accurate expression is (2.7) with the 4th order accurate one given by 



= a. 



■n 




k=l 




(2.14) 



etc. This prescription is easily extended to more general potentials; e.g. 




(2.15) 



is differenced as 



N n-1 




(2.16) 



i=0 k=l 
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For two spatial dimensions, generalization to the Laplacian is trivial. However, potentials 
with cross-terms or different coefficients for d'^ j dx^ and 9^/ dy^ cannot be differenced to 
higher order by this prescription due to the difficulty of eliminating higher derivative terms 
in a multidimensional Taylor expansion. 
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III. The Gowdy Universe Test Case. 

The Gowdy T^ cosmology is conveniently described by the metric [16] 

ds^ = e^/^e^/\- dr^ + dO'^) + e'^ [e^da'^ + 2e^Q dadS + + e"^) db^\ (3.1) 

where A, P, and Q are functions of Q and r only. The T'^ spatial topology is imposed by 
requiring the angular coordinates a and 8 to have arbitrary finite range and < ^ < 27r. 
The time variable r measures the area in the symmetry plane and ^ oo at the singularity 
[9, 12]. The physical interpretation of the polarized model (Q = 0) has been discussed 
extensively [12, 25]. The independent Einstein equations from (3.1) are 

P,,, - - {Q,l - e-2-gJ) = 0, (3.2) 

Q,rr - e-^^Q,ee + 2 (P,. Q,^ - e-^^P,e Q,e) = 0, (3.3) 

X,0 - 2{P,0 P,r + e2^Q,e Q,r ) = 0, (3.4) 

and 

A,. - [P,2 + e-2-P,2 + e2^(Q,2 + e'^^ Q J )] = 0. (3.5) 

The latter two equations for the background X{9,t) are respectively the ^-momentum 
and hamiltonian constraints. [The T"^ topology requires the integral of (3.4) to vanish — 
equivalent to requiring zero total ^-momentum.] 
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Since the "wave" equations (3.2) and (3.3) do not contain A, their evolution is uncon- 
strained. It has been shown, albeit for a different set of variables defined by [16] 

= cosh W + sinh W cos $ 

(3.6) 

e^Q = sinh sin $ 
that the "wave" equations are harmonic map equations for the metric 

dS^ = dP^ + e^^dQ'^. (3.7) 



This is just the harmonic map property of the fields P and Q with (3.7) the metric of the 
target space [15]. The field equations can then be derived from the hamiltonian 

H=lfd9[nl+ e-'^n^Q + e'^^ (Pj + e'^Q,l)] . (3.8) 

Clearly, (3.8) is in the form required by the SI algorithm. 

We note here for future reference that the metric (3.7) admits three Killing fields 
corresponding to the transformations (for constant parameter p) 



P ^P 
Q ^Q + p 



P ^ P-lnp 
Q^pQ 



(3.9) 



(3.10) 



and 



e"" ^ i + Q^) + e-""] + \ [e^(l - Q") - e"^] cosp - e^Qsinp 

e^^Q^e^Qcosp+l [e'P(l-Q2)_ e"^] sinp 
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(3.11) 



This last apparently complicated transformation is just — > W, $ — > $ + ^>o in the 
other coordinates. The presence of the factor e"^"^ in (3.8) and in the wave equations 
suggests that as the singularity at r = oo is approached, the spatial derivatives can be 
neglected yielding the AVTD solution. In the absence of the spatial derivative terms, Eqs. 
(3.2) and (3.3) can be solved exactly in terms of four arbitrary constants a, (3, (, and ^ as 

P = ln[ae-'5^(1 + C^e2^^)] (3.12) 



and 

Substitution of (3.12) and (3.13) into the AVTD form of (3.4) and (3.5) yields 

A = -/3V + Ao. (3.14) 



As r ^ oo, (3.12) and (3.13) become 

P = pT ; Q = Qo (3.15) 



with Qo = l/a( + ^. If a Gowdy solution approaches the AVTD limit, one expects it to 
have the form (3.12)-(3.14) with (in general) different values of a, P, (, and ^ at each 
value of 9. For the polarized case ((5 = 0), it was shown [12] that substitution of (3.15) as 
the limiting form of the exact solution in the metric (3.1) yields the Kasner solution [26] 
with 9 dependent Kasner parameter. In the general case {Q ^ 0), the Kasner axes are 
rotated with respect to the coordinate axes. Isenberg and Moncrief have shown [9] in the 
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polarized case that every solution is AVTD. It is conjectured [13] that this is also true in 
the unpolarized model. 

In the following discussion, we shall consider only the wave equations (3.2) and (3.3) 
since the background X{9,t) may be easily constructed after the dynamical P and Q 
have been found. That P and Q are amplitudes for the two orthogonal polarizations of 
gravitational waves may be seen by analogy with linearized gravity [12, 27]. If the metric 
g^j^L, is expressed as 

9,. = 7i°j + h^^} + (3.16) 
with P and Q assumed small then 

7^ = diag(-eV2g-3r/2^ e^/2e-/2, e-\e--) (3.17a) 

describes a background metric. The designation as background can be enhanced by the 
introduction of spatial averaging [27, 28, 25]. We also find that 

/iW = e-^Ps+^, + e-^ QsxM. (3.176) 
where e+ and Sx are the gravitational wave polarization tensors. In the a-S plane, 

with all other components zero. It is easy to see [12, 27] that in zeroth order, the waves 
act as sources for the background in Eqs. (3.4) and (3.5). In first order, P and Q satisfy 
linear wave equations (3.2) and (3.3). The nonlinearities of the waves enter at the next 
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order. In some sense, we cannot consider these terms to be higher order since the solution 
is quahtatively different if both polarizations are present even if the amplitude of Q is 
small [29]. This is explained by the scaling symmetry (3.10) which implies that a solution 
to (3.2) and (3.3) is independent of the amplitude of Q as long as it is non-zero. 

The polarized case (e.g. Q — 0) yields a single linear wave equation with the general 
exact solution [12] 

oo 

P(d,T) = Zo(ne-^) (an COS nd + bn sin nO) (3.18) 

n=0 

where Zo{x) is a general solution to Bessel's equation of zero order and the a's and 6's are 
arbitrary constants. It is possible to use this exact solution to generate exact "pseudo- 
unpolarized" solutions to the unpolarized equations (3.2) and (3.3) by means of the boost 
symmetry (3.11). Given an exact solution Pq from (3.18), we obtain a class of solutions 

— cosh Pq + sinh Pq cos p 

(3.19) 

e^Q = sinh Pq sin p 

for all values of the parameter < p < 27r. Perhaps the simplest of these solutions is 

P = In cosh Pq 

(3.20) 

Q = tanh Pq 

found for p — tt/2. The pseudo-unpolarized solutions make excellent code tests since 
the fact that they are non- generic is not apparent to the computer. Direct substitution 
of (3.20) into the equations of motion (3.2) and (3.3) shows that nonlinear terms from 
variation of both Hi and H2 [see the discussion in Section II and Eq. (3.8 )] must cancel 
corresponding expressions which arise in the linear terms. Thus the entire code is tested. 
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Grubisic and Moncrief [13] have defined several functions of P and Q and their con- 
jugate momenta which become constant in r in the AVTD regime. They also predict the 
rate of decay to the AVTD regime in terms of the 9 dependent parameters a, and ^. 
For the purposes of this paper, we shall consider only the parameter 

^^(P^+e^^g^)^' (3.21) 

which represents geodesic velocity in the target space (3.7). Grubisic and Moncrief have 
conjectured that for a generic Gowdy T^ model, the AVTD regime will be characterized 
by < V < 1 for all 9. The restriction to generic models must be made because any value 
of V is allowed for polarized solutions [9] and v = 1 is achieved for some non-polarized 
but "asymptotically polarized" solutions [30]. (We note that v is invariant under the 
transformations (3.9)-(3.11) so that pseudo-unpolarized solutions can also have any value 
of V.) The heuristic basis for this conjecture is easily seen. The term e~^'^e^^Q,^ in (3.2) 
becomes e-^^^^-'^^Q,^ in the AVTD limit. Clearly, if v > 1 (with Q,e fixed), this term wiU 
grow, contrary to the AVTD assumption that it is negligible. 
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IV. The Pseudo-Unpolarized Test Case. 

In differenced form, the hamiltonian (3.8) becomes 

i=l 

-2t ^ {A ■\\ 

where (a, 6) = (1/2, 0) and (2/3, —1/24) yields equations correct to second and fourth order 
respectively in the spatial derivatives. [Extension to sixth order requires corresponding 
terms {Pi+3 - PiY, etc. with coefficients (a, 6, c) = (3/4, -3/40, 1/180).] The first sum is 
Hi and the second H2. To evolve with Hi, solve the AVTD solution [(3.12), (3.13), and 
their rJderivatives] for a given P and Q and their conjugate momenta (at each 9 value) for 
the parameters a, /3, and ^. Use these parameters to propagate the initial data to the 
end of the time-step. The evolution with H2 is even easier since it contains no momenta 
so that Pi and Qi remain constant. The momenta are evolved with the (now constant) 
gradients of H2. The overall time dependent factor e"^"^ is then trivially integrated. (We 
note that one may alternatively treat rJ as an extra degree of freedom.) Suzuki's method 
[21] is used to ensure that the time evolution is accurate to the desired order. Greater 
details of our algorithm are given in the Appendix. 

As a code test, initial data appropriate to the pseudo-unpolarized boost (3.20) of the 
exact solution to (3.18) given by 

Po(^,T) = yo(e-")cos^ (4.2) 
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where Yq (x) is an irregular Bessel function were evolved numerically toward the singularity 
[31]. The exact boosted solution is displayed in Fig. 1. Note that the boost transformation 
has generated large 9 derivatives (particularly in Q) due to the hyperbolic tangent. Figures 
2 and 3 illustrate the differences between the numerical and exact solutions for the 2nd 
and 4th order schemes respectively. Although the errors in the second order algorithm are 
small almost everywhere ( 1%), they become unacceptably large (~ 1) as r increases in 
the regions of large 9 derivative. Improvement with the 4th order scheme is dramatic with 
relative errors ~ 10~^ everywhere. 

The accuracy of the 4th order code appears to be acceptable for the following rea- 
sons: The code tests were run with low spatial resolution of 100 total grid points with 99 
(97) representing [0, 27r] for 2nd (4th) order (due to the imposition of periodic boundary 
conditions). The range of r was between and ~ 23. The r interval was chosen for 
convenience — there was no barrier to a much closer approach to the singularity. The large 
spatial gradients (in Q) were correctly represented by the code. We shall see (next section) 
that such features are in fact also characteristic of the generic Gowdy T^ solution. 
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V. Results for a "Generic" Unpolarized Gowdy Model. 

Here we shall discuss the results from a single initial Gowdy T"^ data set for evolution 
toward the singularity. Comparison to other sets appears to indicate that the behavior 
we report here is typical for standing wave initial data. Traveling waves will be discussed 
elsewhere as will the full range of Gowdy phenomenology [18]. Following a suggestion 
by Chrusciel [32], we consider initial data for which the parameter v in (3.20) exceeds 
unity with the initial time dependence that of the AVTD behavior (3.15). The conjecture 
[13] is that a true AVTD regime with v < 1 everywhere will arise during the course of the 
evolution. The selected initial data were 

P = 0, 7rp = 10cos^, Q^cose, ttq = (5.1) 

so that vq = 10 |cos^|. The evolution was performed with 800 spatial grid points for 
the range < r < Gtt with the 4th order accurate code. To study the approach to the 
AVTD regime, we average reported (i.e. saved rather than computed) values of P, Q, and 
V over nearest and next nearest neighbors (to avoid grid scale size structures which must 
be regarded to be unreliable). Figures 4-6 illustrate P, Q, and v respectively vs. 9 for 
various values of r. The approach to the limiting AVTD behavior (3.15) is seen clearly 
in Fig. 7 which shows P, Q, and v vs. r at selected values of 9. The entire evolution is 
displayed in a 3-D surface plot in Figures 8 (for P and Q) and 9 (for v and v < 1). 

We note the following features of the evolution: 
1. The wave amplitude P develops a complicated spatial profile which then freezes after 
which P grows linearly without change of shape. 
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2. The amplitude Q initially grows rapidly [where cos ^ < since there P,,- in (3.3) acts as 
inverse damping] but then becomes constant. (There is some spatial structure in Q which 
does not show up for the amplitude scale used in this graph.) 

3. The parameter v decays essentially monotonically to values < 1 everywhere. This is 
emphasized by Fig. 9b which has a scale adjusted to display only v <1. 

It thus appears that the AVTD regime has been reached and that the parameter v has 
fallen below unity as conjectured. Although we have begun with a single spatial mode, we 
expect the nonlinear terms to cause a generic evolution. Other standing wave initial data 
have been examined. It appears that the evolution is controlled to some extent by fo, the 
initial value of the parameter v. If < 1, there is little growth of spatial structure and the 
AVTD regime is reached quickly. The larger vq, the more nonlinear the wave interactions 
will be. 

If spatial averaging is not used, v > 1 can occur at isolated points where Q,q ~ 0. 
(See Fig. 10.) Such points, in effect, represent the locally polarized models where v > 1 
is allowed. To the extent that Q^e ^ (i.e. that the model really is generic at that 
point), it still evolves to reach the AVTD regime as conjectured at some r >> 67r. The 
spatial averaging dilutes the influence of spiky features in the dynamical variables. The 
development of AVTD behavior with spatial averaging removed is shown in Fig. 11 where 
P/r and v over a limited range in 6 are shown for r = Gtt and r = 147r. The curves (at 
each T value) should be identical in the AVTD limit as r — cxd [see (3.15)]. It appears 
that the asymptotic AVTD behavior will be achieved at sufficiently large r. 

The growth of spatial structure at arbitrarily small scales appears to be characteristic 
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of the Gowdy dynamics for vq >1. Figure 12 reproduces Fig. 8a for < r < 4.2. It is 
easy to see that the initial cos9 spatial profile nonlinearly generates [through (3.2)] cos29 
dependence, etc. The development of this structure ends when the AVTD limit is reached. 
The competition between nonlinear growth and the spatial freezing of AVTD behavior 
suggests that there may exist a vq dependent time-scale to characterize the phenomenon. 
It is possible that this small scale structure may be related to the critical phenomena 
observed by Choptuik [33] and later by Abrahams and Evans [34] for spherical collapse of 
a scalar field and axisymmetric collapse of gravitational waves respectively. 

The presence of this small scale spatial structure which eventually reaches the grid 
spacing scale (unless the AVTD regime is reached) causes the detailed numerical evolution 
to become dependent on the chosen spatial resolution — i.e. at a given r, the finer the grid, 
the smaller the feature that is seen. However, for any r there exists a spatial resolution 
which is sufficient to resolve all the small scale features. This is shown (with no spatial 
averaging) in Fig. 13 which compares the same feature (in this case for ttp) at various 
spatial resolutions. We note that the feature is completely resolved at 6400 grid points 
with a profile that agrees completely with that obtained for 1600 grid points. Greater 
deviations are found for coarser grids. All resolutions represent the solution where it is 
smooth. The evolution of this same feature is shown in Fig. 14. We see that the feature has 
narrowed (and decreased in amplitude) and is no longer resolved at 6400 grid points. This 
narrowing and decreasing amplitude explains the apparent decrease in spatial structure 
seen in the evolution shown in Fig. 4 (although some is due to the increased range in the 
amplitude of P). The spatial averaging used in Figs. 4-9 washes out the structure at small 
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spatial scales. Subsequent evolution of the feature in Fig. 14 shows little change, indicating 
that AVTD behavior (where np ^ constant) is arising. 

This small scale structure is almost certainly a real property of the equations rather 
than a numerical artifact since it can be resolved with sufficient spatial resolution. (As 
a further code test, the structure is seen to disappear when the code is run backward in 
time.) It was also seen in studies of the approach to the singularity using a completely 
different numerical algorithm [35].) Its characterization is currently under investigation. 
Its presence may signal a requirement for adaptive gridding [36] to achieve the spatial 
resolution that appears to be necessary. 
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VI. The Problem. 

It has been shown [14] that an arbitrary cosmological spacetime onT^ x R containing 
a spacehke ^7(1) symmetry can be described by the 5 degrees of freedom (fi,u;,x,z,A and 
their respective conjugate momenta p, r, Px,Pz,PA- All variables are functions of the spatial 
coordinates u and v and a time coordinate r which measures the size of the universe in the 
symmetry direction. The conformal 2-metric of the space orthogonal to the Killing field is 

^ /e^'+ e-2^(l + x)^ e^^ + e'^'ix^ - 1) \ 
eab = ^ (6.1) 

with a,h — 1,2 and det Cab — 1- (The 2-metric itself is Qah — ^ah-) The scalar curvature 
of this conformal 2-space is 



X, 



. 2z , —22 I —2z _ —2z 2 ] 



(6.2) 



+ [-6-^^(1 + x)x,y +e-^^(l + xfz,y -e^^z^y +e-^^x, 



u 



I 2z I —22 _|_ —22 _ —22 2 ] 



Given these definitions, the field equations (for the dynamical variables in the gauge N = 
= where N is the lapse and '^g is the determinant of the 2-metric) can be derived 
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from the hamiltonian 

H = - j jdudv {\pI + \e^^pl + 1/ + \e^^r^ - \pl - 2p^) 

j du dv [-e^i? + (e«''A,a) ,6 +2e^e«V,a ¥^,6 +|e^e-4'^e«^a;,a a;,^] . 

(6.3) 

We note that the integrand in Eq. (6.3) is not the superhamiltonian, but differs from it by 
an overall sign and an additional term linear in p\. The sign arises from the fact that r 
results from an original time variable t via dt = — e~'^dr while the additional term comes 
from the fact that a time dependent transformation is required to obtain our variables 
from the original Arnowitt, Deser, and Misner (ADM) [37] ones. 

It is clear that the hamiltonian (6.3) fits naturally into the form of the SI operator 
splitting. Even more striking is the fact that the first integral, Hi, contains two copies of 
the Gowdy kinetic term [see (3.8)] plus the kinetic energy of a free particle (with the 
"wrong" sign). Thus we already have the exact solution for Hi from (3.12) and (3.13) as 
in the Gowdy case plus the trivial free particle solution. The second integral, H2i is more 
difficult — not to solve since there are no momenta so all the variables are assumed constant 
over the sub-time step — but to spatially difference to the correct order of accuracy. 

In addition, the dynamical degrees of freedom are constrained — the Gowdy split into 
wave and background variables does not occur. For example, the integrand of (6.3), leaving 
out the term linear in p^, is proportional to the hamiltonian constraint. There are also non- 
trivial momentum constraints. Since the same terms occur in H2 and in the hamiltonian 
constraint, it is probably advantageous to preserve the divergence structure of the first 
two terms in the square bracket rather than to partially integrate them. A differencing 
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scheme for the field equations that recognizes the underlying structure of the constraints 
may aid in keeping the numerical evolution on the constraint hypersurface. Although the 
differencing scheme for H2 outlined in Section II can be extended to the Laplacian, it 
becomes problematical for the mixed partial derivatives that will arise here. We plan to 
begin with a plausible 2nd order accurate scheme. 

Since the dynamical evolution is constrained, specification of initial data is non-trivial. 
Fortunately, examples of "half-polarized" exact solutions of the constraint equations are 
known [38]. 

Although almost nothing is known about the generic behavior of the U{1) models, 
it is known that not all solutions can be AVTD. This is because this class of solutions 
contains the Mixmaster cosmologies [2- 4] which (since the influence of the potential never 
disappears) are not AVTD. If, as has been conjectured, the Mixmaster dynamics cannot 
survive the presence of the spatial inhomogeneity degrees of freedom, the models could 
still be AVTD. We note here, however, that the transformation necessary to obtain our 
variables has obscured the meaning of AVTD since the transformation to the twist potential 
degree of freedom [a;, r] has interchanged coordinates and momenta. It is also expected 
that small scale spatial structure will appear in the generic case. The possible need for high 
spatial resolution in two spatial dimensions may push the limits of current computational 
technology. 
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VII. Summary and Conclusions. 

A SI scheme that sphts the hamiltonian into exactly solvable kinetic and potential 
pieces is ideally suited to the numerical study of the singularity structure of spatially in- 
homogeneous cosmologies. For both the Gowdy model and the more general U{1) 
problem (on x R), the dynamical equations arise from a variational principle that also 
splits naturally into kinetic and potential pieces which separately can be solved exactly. 
The exact solution for the kinetic sub-hamiltonian is in fact just the AVTD solution con- 
jectured to arise for the generic Gowdy T'^ model (taken twice plus a free particle term in 
the U{1) problem). 

The SI code has been implemented through 4th order in both time (using Suzuki's 
method) and space for the generic Gowdy model. Comparison for a pseudo-unpolarized 
test case (obtained by a boost symmetry from an exact solution for the polarized model) 
shows agreement to 1 : 10^ for the 4th order accurate code. The 2nd order accurate code 
diverges unacceptably in regions of large spatial derivative as the singularity is approached. 
The cost of the extra accuracy is essentially a factor of three in computational time (since 
three 2nd order time steps are required to produce a 4th order time step). The extra 
spatial accuracy involves negligible cost. 

For reasonably generic Gowdy initial data, we have been able to support the conjecture 
[13] that the models are AVTD with v < 1 everywhere. More detailed study of the approach 
to the AVTD regime to compare with the detailed predictions in [13] is in progress. 

An interesting new phenomenon of the development of small scale spatial structure has 
been observed. Studies to characterize this behavior in terms of the competition between 
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nonlinear generation of short wavelength modes and the freezing of the spatial profile in 
the AVTD regime are underway. 

Thus we have applied the SI scheme to the unpolarized Gowdy cosmology and have 
been able to test the code, to study and verify AVTD regime conjectures, and to discover 
the new phenomenon of nonlinear small scale structure. An even richer phenomenology 
awaits the application of this method to the unknown territory of the U(l) problem. 
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Appendix: Details of the Second Order Accurate SI Algorithm for the Gowdy Model. 
The hamiltonian (3.8) can be put in differenced form 



TV 



N 



Q 



,) 



i=0 



~ I Pi+i+Pi ( Qi+'^ — Q.^ 



{A.l) 



with each sum regarded to be an independent sub-hamiltonian [Hi and H2 respectively). 
Given initial data Qi{T^), Pi{T^),7^Q^{T^),7^p^{T^), we use H2 to evolve to + V2AT to 
yield 

Qi = Qiir'), 



Pi = Pi{r'), 



■2r,- 



dV 



dQi 



{A.2) 



7fp, = TTp,{t^) + 36 



1 „-2r,- 



-(-^j+i-Tj) _ I 



dV 



J dP, 



where V is the "spatially dependent" part of the potential term in (A.l). Note that the 
time dependence has been taken into account separately and that two terms in the sum 
contribute to the indicated gradient. 

Now solve the AVTD equations (3.12) and (3.13) and their r derivatives using Qi, Pi, 
ttq. , ^tp^ for the constants Q, Pi, ai, and ^i . We find 
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1. If TTp^ ^ and ttq^ ^ then 



1 + C 

Cf-1 



p. 

e " 



(A.3a) 



2. If 7fp. ^ 6wt Ti-Q, = then 

Ci = 0, /3i = -Tfp,, ai = e-^S = Qi. (^.4) 

3. If TTp. =0 then for any ttq^, 

Ci = 1, /3i = -^Q,e-^S 



(A.5) 



These values are then used to evolve the Qi, -Pi, ttq., 7fp. with the AVTD equations to 
+ At. This is just the evolution by Hi required by the algorithm. We find 

7rQi=7rQ,, {A.6a) 
29 



{AM) 



di e 



-/3i(Tj+i-rj) 



_|_ ^2g/3i(rj + i-rj)^ 



(A.6c) 



„2/3i(rj + i-Tj) 



(yl.6d) 



Finally, the Qi, Pi,'7^Q^,7tp^ are evolved with for At/2 to yield the original variables 
updated to the next time-step: 



dV 



7rp,{T^+^) = np, + ie 



2^ 



ay 



Qk,Pk 



Qk,Pk 



{A.7) 



To achieve fourth order accuracy in time, repeat the entire procedure three times with 
steps of sAt, (1 — 2s)At, and sAr respectively with s = (2 — 2-^/^)"^. For fourth order 
spatial accuracy, replace the potential term in (A.l) by the appropriate version of (4.1). 
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FIGURE CAPTIONS 



Fig. 1. The exact solution for the pseudo-unpolarized test case. Initial data are generated 
from Eq. (3.20) applied to (4.2) (and its r derivative) evaluated at r = 0. The axis scales 
for 9 and r are [0, 2tt\ and [0, 23] respectively. The vertical axis indicates the values of (a) 
P(6',t) and (b) g(6',T). 

Fig. 2. Errors in the 9-t plane for the 2nd order accurate numerical scheme. The ranges 
of 9 and T are the same as in Fig. 1. (a) The vertical scale measures Pnumencal — Pexact 
(where Pexact is shown in Fig. 1(a)) and Pnumericai is the computation, (b) The same as 
(a) but for Q. 

Fig. 3. Errors in the 9-t plane for the 4th order accurate numerical scheme. All axes and 
variables are the same as in Fig. 2 except that the 4th order accurate numerical scheme 
has been used. 

Fig. 4. P{9,t) vs. 9 at selected values of r. Results of the simulation for initial data 
from (5.1) in the 4th order accurate numerical scheme are shown. The computation was 
performed with 800 spatial grid points. The display consists of values at 100 of these grid 
points averaged over nearest and next nearest neighbors in the full array. (The selection 
of spatial grid points accounts for any asymmetry about 9 = tt.) Six graphs at selected r 
values are stacked. In all cases, the horizontal axis is < ^ < 27r. The numerical scales 
on either the left or right axis denote the amplitude of P. The simulation represents the 
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range < r < 18.85. 

Fig. 5. QJ{9, r) vs. 9 at selected values of r. The same as Fig. 4 but for Q. For r > 1.14, 

the vertical scale is [0, 500]. 

Fig. 6. v{9,t) vs. 9 at selected values of r. The same as Fig. 4 but for v as defined by 
(3.21). 

Fig. 7. P{9^ r), Q{9^ r), and v{9^ r) vs. r at selected values of 9. In all cases, the horizontal 
axis is < T < 18. On the graph for Q, the values of the solid lines correspond to the left 
axis and the dashed lines to the right axis. 

Fig. 8. P and Q in the 9-t plane. The complete results of the simulation in Figs. 4-7 are 
shown in the 9-t plane. The axis scales for 9 and r are [0, 2tt] and [0, 67r] respectively. The 
vertical axes denote the values of (a) P and (b) Q respectively. 

Fig. 9. The values of v in the 9-t plane. This figure is the same as Fig. 8 but for the 
parameter v. Note that the viewing angle has changed for greater clarity, (a) All values 
of V . (b) The same data with the vertical scale chosen to display only t> < 1. 

Fig. 10. A "non-generic" point in 9. Q and f vs. 6* at r = 6n are shown near a point 
with V > 1 for the initial data of Figs. 4-9 for a simulation with 6400 spatial grid points 
with no spatial averaging of the results. The graphs are produced using all the spatial 
resolution available in the simulation. The left axis corresponds to v and the right to Q. 

36 



Note that v > 1 only where Q,g ^ and that the ranges displayed for Q and 9 are small. 
The horizontal dashed line denotes v = 1. 

Fig. 11. Approaching the AVTD limit. Graphs of v (solid line) and P/r (dashed line) 
vs. 9 are overlaid. The data are taken from the same simulation as Figs. 4-9 for 800 spatial 
grid points but with no spatial averaging. The vertical scale for both P/t and v ranges 
from to 1.2 while the horizontal axis is 2.6 < ^ < 3.2. (a) r = 18.86; (b) r = 43.96. 

Fig. 12. Growth of small scale spatial structure. The graph shows the data for P{9, r) 
from Fig. 8(a) for the range < r < 4.2. 

Fig. 13. Spatial resolution dependence of generic spiky features. A typical spiky feature in 
TTp vs. ^ at T = 27r is shown at resolutions of 400, 800, 1600, and 6400 spatial grid points. 

Fig. 14. Evolution of a spiky feature. The dashed line indicates the same plot as Fig. 13 at 
T = 27r for 6400 spatial grid points on the expanded scale of this figure. The same feature 
at 6400 spatial grid points is shown at r = 47r (solid line) and r = 67r (dashed line with 
circles). For comparison, the feature computed with 1600 spatial grid points is shown at 
r = Att (dashed line with squares). 
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